The coral growth data for this analysis variables related to the growth of individual colonies collected from the visual analysis of x-rays of 1 inch wide slices cut down the central growth plane of individual colonies from three Persian Gulf sites and two sites off the coast of Oman. The growth variables are measured in yearly increments, so environmental data must also be summarized in yearly increments.
Lets look at the structure of the growth data first:
onyr <- read.csv("rawdata/maxvalues-year2000onwards.csv")
head(onyr)
## year site colony track sample calc extension density hdbands partmort
## 1 2000 AlAqah ala12 t1 ala12_t1 0.797 0.6559 1.534612 0 0
## 2 2000 AlAqah ala12 t2 ala12_t2 1.007 0.5701 1.397308 0 0
## 3 2000 Dibba dib01 t1 dib01_t1 0.956 0.7467 1.280707 0 0
## 4 2000 Dibba dib01 t2 dib01_t2 0.526 0.3784 1.390043 0 0
## 5 2000 Dibba dib04 t1 dib04_t1 0.940 0.8173 1.150707 0 0
## 6 2000 Dibba dib04 t2 dib04_t2 1.335 1.0998 1.213488 0 0
## height dia
## 1 13.5 23
## 2 13.5 23
## 3 NA NA
## 4 NA NA
## 5 14.0 28
## 6 14.0 28
summary(onyr)
## year site colony track
## Min. :2000 Length:768 Length:768 Length:768
## 1st Qu.:2005 Class :character Class :character Class :character
## Median :2008 Mode :character Mode :character Mode :character
## Mean :2008
## 3rd Qu.:2011
## Max. :2013
##
## sample calc extension density
## Length:768 Min. :0.1650 Min. :0.1665 Min. :0.8693
## Class :character 1st Qu.:0.5965 1st Qu.:0.4440 1st Qu.:1.1948
## Mode :character Median :0.7695 Median :0.5750 Median :1.3423
## Mean :0.8211 Mean :0.6044 Mean :1.3722
## 3rd Qu.:1.0170 3rd Qu.:0.7232 3rd Qu.:1.5187
## Max. :1.9930 Max. :2.9767 Max. :2.3394
## NA's :4 NA's :9 NA's :7
## hdbands partmort height dia
## Min. :0.00000 Min. :0.000000 Min. : 4.000 Min. :11.00
## 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.: 7.000 1st Qu.:15.00
## Median :0.00000 Median :0.000000 Median : 8.500 Median :16.00
## Mean :0.08333 Mean :0.003906 Mean : 8.814 Mean :17.94
## 3rd Qu.:0.00000 3rd Qu.:0.000000 3rd Qu.:11.000 3rd Qu.:22.00
## Max. :1.00000 Max. :1.000000 Max. :14.000 Max. :28.00
## NA's :50 NA's :50
# number of unique colonies
length(unique(onyr$colony))
## [1] 38
# number of unique sites = 5
unique(onyr$site)
## [1] "AlAqah" "Dibba" "RasGhanada" "Saadiyat" "Delma"
# plot of individual colony diameter Vs height
xyplot(height ~ dia | site, data = onyr,
#group = year, #t="l",
grid = TRUE, auto.key=list(space="top", columns=5,
title="Site", cex.title=1,
lines=TRUE, points=FALSE))
So we have 768 measurements from 38 colonies from five sites. The coral growth data record growth from 2000-2013.
Lets look at the structure of the environmental variables, which are data products derived from satellite observations. They consist of sea surface temperature metrics and current speed and direction metrics from each of the sites.
Here is the Sea Surface Temperature (SST) data for each site. We can see that the data product summarizes the SST daily and the dataset starts in 1985:
temp <- read.csv("rawdata/ctemp_sst_mod.csv")
str(temp)
## 'data.frame': 12370 obs. of 11 variables:
## $ date : chr "1985JAN01" "1985JAN02" "1985JAN03" "1985JAN04" ...
## $ date1 : chr "1/1/85" "2/1/85" "3/1/85" "4/1/85" ...
## $ yearday : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Year : int 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
## $ Delma : num 22.7 22.6 22.5 22.2 22.1 ...
## $ RasGhanada: num 23.3 23.1 23 22.8 22.8 ...
## $ Saadiyat : num 23.3 23.1 23 22.7 22.7 ...
## $ AlAqah : num 24 23.9 23.9 23.9 23.9 ...
## $ Dibba : num 24 23.9 23.9 23.8 23.8 ...
## $ month : chr "JAN" "JAN" "JAN" "JAN" ...
## $ season : chr "winter" "winter" "winter" "winter" ...
## Resetting variable names to make it easier to reshape dataset after adding moving averages
temp1 <- temp[,5:9]
colnames(temp1) <- paste(colnames(temp1), "raw", sep = ".")
temp <- cbind(temp[,c(1,2,3,4,10,11)], temp1)
Now lets smooth out some the the effect of daily variation and get some moving averages of temporal trends. We will use a centred moving average to capture coldest and hottest periods in each year:
library(tidyverse)
temp11 <- temp %>%
select(date, date1, yearday, Year, month, season, Delma.raw, RasGhanada.raw, Saadiyat.raw, AlAqah.raw, Dibba.raw) %>%
mutate(Delma.CMA = rollmean(Delma.raw, k = 61, fill = NA),
RasGhanada.CMA = rollmean(RasGhanada.raw, k = 61, fill = NA),
Saadiyat.CMA = rollmean(Saadiyat.raw, k = 61, fill = NA),
AlAqah.CMA = rollmean(AlAqah.raw, k = 61, fill = NA),
Dibba.CMA = rollmean(Dibba.raw, k = 61, fill = NA))
## Reformatting and making some plots to visualize SST. Only looking at SST from one site and two different years:
##
temp11$date1 <- as.Date(temp11$date1, '%d/%m/%y')
str(temp)
## 'data.frame': 12370 obs. of 11 variables:
## $ date : chr "1985JAN01" "1985JAN02" "1985JAN03" "1985JAN04" ...
## $ date1 : chr "1/1/85" "2/1/85" "3/1/85" "4/1/85" ...
## $ yearday : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Year : int 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
## $ month : chr "JAN" "JAN" "JAN" "JAN" ...
## $ season : chr "winter" "winter" "winter" "winter" ...
## $ Delma.raw : num 22.7 22.6 22.5 22.2 22.1 ...
## $ RasGhanada.raw: num 23.3 23.1 23 22.8 22.8 ...
## $ Saadiyat.raw : num 23.3 23.1 23 22.7 22.7 ...
## $ AlAqah.raw : num 24 23.9 23.9 23.9 23.9 ...
## $ Dibba.raw : num 24 23.9 23.9 23.8 23.8 ...
require(ggplot2)
temp98 <- subset(temp11, Year ==2002 | Year ==2003)
ggplot( data = temp98, aes(date1, Delma.raw )) + geom_point() +
#geom_line(data = temp98, aes(date1, Delma.raw ), color = "red")+
geom_line(data = temp98, aes(date1, Delma.CMA ), color = "blue")
The data needs reshaping so raw SST (Delma, RasGhanada, Saadiyat, AlAqah, Dibba) and 61 day centered averaged temps are in single columns with a site factor level instead (removing first 30 & last 30 days of data first to remove all NAs produced from the centred moving average)
temp11 <- temp11[31:(nrow(temp11)-31),]
templong <- temp11 %>%
gather("key", "value", -date, -date1, -yearday, -Year, -month, -season) %>%
separate(key, c("Site", "raw_CMA"), sep = "\\.") %>%
spread(raw_CMA, value)
xyplot(raw ~ yearday | as.factor(Year) , data = templong,
group = Site, t="l",
grid = TRUE, auto.key=list(space="top", columns=5,
title="Site", cex.title=1,
lines=TRUE, points=FALSE))
So we can see the expected annually cyclical nature of SST here.
Because we are interested in summarizing variables for each calendar year, we need a way to discriminate the different seasons from the temperature data itself in order to summarize, say seasonal(summer / winter) average temperatures or to use the rates of change between summer highs to winter and winter lows to summer each year. We need to do this because using static calendar dates to define these may not capture the appropriate data range for particular summary, because of temporal variation in seasons (e,g, a “late summer” or “early winter”. To begine with, lets detect and plot the90 days after coldest (blue) and hottest day of each year are detected (red) - checking a few different sites & only looking for a few different years 2003-2008
plot_df <- templong %>%
mutate(year = lubridate::year(date1)) %>%
group_by(year, Site) %>%
mutate(raw_min = +(raw == min(raw)),
raw_max = +(raw == max(raw))) %>%
ungroup() %>%
mutate(raw_min = cumsum(raw_min - lag(raw_min, 90, default = 0)),
raw_max = cumsum(raw_max - lag(raw_max, 90, default = 0)))
ggplot(
plot_df %>%
filter(Site == "Delma", Year >= 2003, Year <= 2008),
aes(date1, raw)
) +
geom_line() +
geom_vline(
aes(xintercept = date1),
data = plot_df %>% filter(Site == "Delma", Year >= 2003, Year <= 2008, raw_min > 0),
alpha = 0.1,
colour = "blue"
) +
geom_vline(
aes(xintercept = date1),
data = plot_df %>% filter(Site == "Delma", Year >= 2003, Year <= 2008, raw_max > 0),
alpha = 0.1,
colour = "red"
)
ggplot(
plot_df %>%
filter(Site == "Dibba", Year >= 2003, Year <= 2008),
aes(date1, raw)
) +
geom_line() +
geom_vline(
aes(xintercept = date1),
data = plot_df %>% filter(Site == "Dibba", Year >= 2003, Year <= 2008, raw_min > 0),
alpha = 0.1,
colour = "blue"
) +
geom_vline(
aes(xintercept = date1),
data = plot_df %>% filter(Site == "Dibba", Year >= 2003, Year <= 2008, raw_max > 0),
alpha = 0.1,
colour = "red"
)
We can now subset the raw data to be from the yearly minimum SST at each site for the 90 days after.
raw_minroc <- subset(plot_df, raw_min > 0)
xyplot(raw ~ yearday | as.factor(Year) , data = raw_minroc,
group = Site, type ="l", xlim = c(0,200),
grid = TRUE, auto.key=list(space="top", columns=5,
title="Site", cex.title=1,
lines=TRUE, points=FALSE))
Ideally if we fit a simple linear regression to this line, it would give a measure of the rate of change toward summer from the coldest day of winter in that year. But it seems like a few problems to deal with first - For example, 2005 something weird has gone on at the site RasG. Other instances is where the coldestr day is on the 31st December, (RasG and Saadiyat) which bugs up the by year filtering, so we need to remove (RasG and Saadiyat, 31 dec 2004) running filtering again without those rows
templong <- templong %>%
filter(date1 != "2004-12-31" | Site != "RasGhanada" ) %>%
filter(date1 != "2004-12-31" | Site != "Saadiyat")
nrow(templong)
## [1] 61543
plot_df <- templong %>%
mutate(year = lubridate::year(date1)) %>%
group_by(year, Site) %>%
mutate(raw_min = +(raw == min(raw)),
raw_max = +(raw == max(raw))) %>%
ungroup() %>%
mutate(raw_min = cumsum(raw_min - lag(raw_min, 90, default = 0)),
raw_max = cumsum(raw_max - lag(raw_max, 90, default = 0)))
str(plot_df)
## tibble [61,543 × 12] (S3: tbl_df/tbl/data.frame)
## $ date : chr [1:61543] "1985JAN31" "1985FEB01" "1985FEB02" "1985FEB03" ...
## $ date1 : Date[1:61543], format: "1985-01-31" "1985-02-01" ...
## $ yearday: int [1:61543] 31 32 33 34 35 36 37 38 39 40 ...
## $ Year : int [1:61543] 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
## $ month : chr [1:61543] "JAN" "FEB" "FEB" "FEB" ...
## $ season : chr [1:61543] "winter" "winter" "winter" "winter" ...
## $ Site : chr [1:61543] "Delma" "Delma" "Delma" "Delma" ...
## $ CMA : num [1:61543] 20.9 20.9 20.8 20.8 20.7 ...
## $ raw : num [1:61543] 21.1 20.9 21 21 21.2 ...
## $ year : num [1:61543] 1985 1985 1985 1985 1985 ...
## $ raw_min: int [1:61543] 0 0 0 0 0 0 0 0 0 0 ...
## $ raw_max: int [1:61543] 0 0 0 0 0 0 0 0 0 0 ...
Now lets plotting again to make sure 90 days after coldest (blue) and hottest day of each year are detected (red) - checking a few different sites
ggplot(plot_df %>% filter(Site =="Delma"), aes(date1, raw)) +
geom_line() +
geom_vline(aes(xintercept = date1), plot_df %>% filter(Site =="Delma", raw_min > 0),
alpha = 0.1, colour = "blue") +
geom_vline(aes(xintercept = date1), plot_df %>% filter(Site =="Delma", raw_max > 0),
alpha = 0.1, colour = "red")
ggplot(plot_df %>% filter(Site =="Dibba"), aes(date1, raw)) +
geom_line() +
geom_vline(aes(xintercept = date1), plot_df %>% filter(Site =="Dibba", raw_min > 0),
alpha = 0.1, colour = "blue") +
geom_vline(aes(xintercept = date1), plot_df %>% filter(Site =="Dibba", raw_max > 0),
alpha = 0.1, colour = "red")
Now lets check to see if the data are now amenable to generating yearly rates of change between seasons for the 90 days after the coldest and hottest day of each year.
##checking
raw_minroc <- subset(plot_df, raw_min > 0)
raw_maxroc <- subset(plot_df, raw_max > 0)
## winter to summer
xyplot(raw ~ yearday | as.factor(Year) , data = raw_minroc,
group = Site, type ="l", xlim = c(0,200),
grid = TRUE, auto.key=list(space="top", columns=5,
title="Site", cex.title=1,
lines=TRUE, points=FALSE))
##Summer to winter
xyplot(raw ~ yearday | as.factor(Year) , data = raw_maxroc,
group = Site, type ="l", #xlim = c(0,200),
grid = TRUE, auto.key=list(space="top", columns=5,
title="Site", cex.title=1,
lines=TRUE, points=FALSE))
There appears to be a problem related to when the data starts for sites RasG and Al Aqah for summer to winter in 1985 but this will be discarded anyway (since coral data is only from 2000).
In order to gather a measure of seasonal rates of change, we will fit a linear model within each year for SST~yearday where yearday will be normalized to the coldest (or hottest) day to day 1 before fitting regression so that intercepts are normalised to hottest/coldest day (as applicable)
library(data.table)
str(raw_minroc)
## tibble [15,304 × 12] (S3: tbl_df/tbl/data.frame)
## $ date : chr [1:15304] "1985MAR08" "1985MAR09" "1985MAR10" "1985MAR11" ...
## $ date1 : Date[1:15304], format: "1985-03-08" "1985-03-09" ...
## $ yearday: int [1:15304] 67 68 69 70 71 72 73 74 75 76 ...
## $ Year : int [1:15304] 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
## $ month : chr [1:15304] "MAR" "MAR" "MAR" "MAR" ...
## $ season : chr [1:15304] "spring" "spring" "spring" "spring" ...
## $ Site : chr [1:15304] "Delma" "Delma" "Delma" "Delma" ...
## $ CMA : num [1:15304] 20.6 20.6 20.6 20.7 20.7 ...
## $ raw : num [1:15304] 18.8 18.8 18.9 19 18.9 ...
## $ year : num [1:15304] 1985 1985 1985 1985 1985 ...
## $ raw_min: int [1:15304] 1 1 1 1 1 1 1 1 1 1 ...
## $ raw_max: int [1:15304] 0 0 0 0 0 0 0 0 0 0 ...
str(raw_maxroc)
## tibble [15,281 × 12] (S3: tbl_df/tbl/data.frame)
## $ date : chr [1:15281] "1985SEP13" "1985SEP14" "1985SEP15" "1985SEP16" ...
## $ date1 : Date[1:15281], format: "1985-09-13" "1985-09-14" ...
## $ yearday: int [1:15281] 256 257 258 259 260 261 262 263 264 265 ...
## $ Year : int [1:15281] 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
## $ month : chr [1:15281] "SEP" "SEP" "SEP" "SEP" ...
## $ season : chr [1:15281] "fall" "fall" "fall" "fall" ...
## $ Site : chr [1:15281] "Delma" "Delma" "Delma" "Delma" ...
## $ CMA : num [1:15281] 33.3 33.3 33.2 33.2 33.2 ...
## $ raw : num [1:15281] 35.7 35 34.4 34.6 34.4 ...
## $ year : num [1:15281] 1985 1985 1985 1985 1985 ...
## $ raw_min: int [1:15281] 0 0 0 0 0 0 0 0 0 0 ...
## $ raw_max: int [1:15281] 1 1 1 1 1 1 1 1 1 1 ...
raw_minroc$Site <- as.factor(raw_minroc$Site)
raw_maxroc$Site <- as.factor(raw_maxroc$Site)
raw_minroc <- raw_minroc %>%
group_by(Site, year) %>%
arrange(date1) %>%
mutate(normday = yearday - min(yearday) + 1)
raw_maxroc <- raw_maxroc %>%
group_by(Site, year) %>%
arrange(date1) %>%
mutate(normday = yearday - min(yearday) + 1)
##
library(dplyr)
library(tidyr)
library(purrr)
library(broom)
# Fit winter models by Site and year for SST rates of change:
wintrocfin <- raw_minroc %>%
group_by(Site, year) %>%
nest() %>%
mutate(
model = map(data, ~ lm(raw ~ normday, data = .x)),
tidied = map(model, tidy),
glanced = map(model, glance)
) %>%
unnest(cols = c(tidied), names_sep = "_") %>%
filter(tidied_term == "normday") %>%
mutate(
winSlope = tidied_estimate,
winIntercept = map_dbl(model, ~ coef(.x)[1]),
winRSquared = map_dbl(glanced, "r.squared")
) %>%
select(Site, year, winRSquared, winIntercept, winSlope) %>%
left_join(
raw_minroc %>%
group_by(Site, year) %>%
summarise(ccoldday = min(yearday), .groups = "drop"),
by = c("Site", "year")
)
# Fit summer models by Site and year for SST rates of change:
summrocfin <- raw_maxroc %>%
group_by(Site, year) %>%
nest() %>%
mutate(
model = map(data, ~ lm(raw ~ normday, data = .x)),
tidied = map(model, tidy),
glanced = map(model, glance)
) %>%
unnest(cols = c(tidied), names_sep = "_") %>%
filter(tidied_term == "normday") %>%
mutate(
sumSlope = tidied_estimate,
sumIntercept = map_dbl(model, ~ coef(.x)[1]),
sumRSquared = map_dbl(glanced, "r.squared")
) %>%
select(Site, year, sumRSquared, sumIntercept, sumSlope) %>%
left_join(
raw_maxroc %>%
group_by(Site, year) %>%
summarise(warmday = max(yearday), .groups = "drop"),
by = c("Site", "year")
)
pairs(summrocfin[,3:6])
Now that we have temperature rates of changed summarized, lets move onto to summarizing the hottest and coldest 61 days of each year to indicate the relative “hotness” of the annual summer and “coldness” of the annual winter at each site. For this we will use moving averages rather than daily temperatures, so as to remove daily noise from finding the hottest and coldest periods of each year.
library(dplyr)
library(lubridate)
library(purrr)
# Adjust year so that DECEMBER is counted as part of the following winter - this is for years where the coldest period is close to the start of the year and therefore an attempt to extract the previous 30 days from the coldest moving average value in each yea is not possible.
templong <- templong %>%
mutate(yearAD = ifelse(month == "DEC", Year + 1, Year))
# Identify coldest (min CMA) and hottest (max CMA) days per site-year
sumwin_df <- templong %>%
group_by(yearAD, Site) %>%
mutate(
CMA_min = as.integer(CMA == min(CMA, na.rm = TRUE)),
CMA_max = as.integer(CMA == max(CMA, na.rm = TRUE))
) %>%
ungroup()
# Helper function to extract 61-day window around target indices
extract_window <- function(df, flag_col) {
inds <- which(df[[flag_col]] == 1)
rows <- unique(unlist(lapply(inds, function(i) (i - 30):(i + 30))))
rows <- rows[rows > 0 & rows <= nrow(df)] # Keep only valid indices
result <- df[rows, ]
result[[flag_col]] <- 1 # Mark entire window with flag
return(result)
}
# Extract 61-day windows for coldest (winter) and hottest (summer) days
sumwin_win <- extract_window(sumwin_df, "CMA_min") %>%
mutate(yearAD = ifelse(month == "DEC", Year + 1, Year)) # Ensure DEC counted in correct year
sumwin_sum <- extract_window(sumwin_df, "CMA_max") # Summer doesn't need year fix
# Merge window flags back into full dataset
temp1 <- templong %>%
left_join(sumwin_sum %>% select(date1, Site, CMA_max), by = c("date1", "Site")) %>%
left_join(sumwin_win %>% select(date1, Site, CMA_min), by = c("date1", "Site"))
Lets now have a look at a few sites and time period to ensure we are indeed extracting the hottest and coldest periods of each season in each year at each site
We can also visualise varition between sites
Looks good. Now we can easily summarize to extract minimum, maximum, range, mean and the SD (standard deviation) in temperature (representing temperature variability) for the hottest and coldest 61 days of each year at each site.
library(dplyr)
library(purrr)
# Create winteryear column for situations where the coldest period overlaps with the start of the year
temp1 <- temp1 %>%
mutate(
Site = as.factor(Site),
winteryear = ifelse(month == "DEC", Year + 1, Year)
)
# Subset winter and summer temperature windows
wtemp <- temp1 %>% filter(CMA_min == 1)
stemp <- temp1 %>% filter(CMA_max == 1)
# Generic summarisation function for SST data
summarise_sst <- function(data, group_var, prefix) {
data %>%
group_by(across(all_of(group_var)), Site) %>%
summarise(
maxM = max(CMA, na.rm = TRUE),
minM = min(CMA, na.rm = TRUE),
avM = mean(CMA, na.rm = TRUE),
sdM = sd(CMA, na.rm = TRUE),
maxR = max(raw, na.rm = TRUE),
minR = min(raw, na.rm = TRUE),
avR = mean(raw, na.rm = TRUE),
sdR = sd(raw, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(site = Site) %>%
select(-Site) %>%
rename_with(~ paste0(prefix, .), -c(all_of(group_var), site))
}
# Get summaries
Wtsumall <- summarise_sst(wtemp, "winteryear", "Wt") %>%
rename(year = winteryear) %>%
mutate(
WtrangeM = WtmaxM - WtminM,
WtrangeR = WtmaxR - WtminR
) %>%
drop_na()
Stsumall <- summarise_sst(stemp, "Year", "St") %>%
rename(year = Year) %>%
mutate(
StrangeM = StmaxM - StminM,
StrangeR = StmaxR - StminR
) %>%
drop_na()
Finally lets combine these temperature summaries with the coral growth data:
# Prepare coral growth dataset
onyr <- onyr %>%
mutate(
site = as.factor(site),
Site = site
)
# Rename 'site' columns in SST rate-of-change summaries
wintrocfin <- wintrocfin %>%
rename(site = Site)
summrocfin <- summrocfin %>%
rename(site = Site)
# Merge everything
test <- onyr %>%
left_join(Wtsumall, by = c("year", "site")) %>%
left_join(Stsumall, by = c("year", "site")) %>%
left_join(wintrocfin, by = c("year", "site")) %>%
left_join(summrocfin, by = c("year", "site"))
str(test)
## 'data.frame': 768 obs. of 41 variables:
## $ year : num 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
## $ site : Factor w/ 5 levels "AlAqah","Delma",..: 1 1 3 3 3 3 4 4 4 4 ...
## $ colony : chr "ala12" "ala12" "dib01" "dib01" ...
## $ track : chr "t1" "t2" "t1" "t2" ...
## $ sample : chr "ala12_t1" "ala12_t2" "dib01_t1" "dib01_t2" ...
## $ calc : num 0.797 1.007 0.956 0.526 0.94 ...
## $ extension : num 0.656 0.57 0.747 0.378 0.817 ...
## $ density : num 1.53 1.4 1.28 1.39 1.15 ...
## $ hdbands : int 0 0 0 0 0 0 0 0 0 0 ...
## $ partmort : int 0 0 0 0 0 0 0 0 0 0 ...
## $ height : num 13.5 13.5 NA NA 14 14 7.5 7.5 12 12 ...
## $ dia : int 23 23 NA NA 28 28 14 14 16 16 ...
## $ Site : Factor w/ 5 levels "AlAqah","Delma",..: 1 1 3 3 3 3 4 4 4 4 ...
## $ WtmaxM : num 24.1 24.1 24 24 24 ...
## $ WtminM : num 23.1 23.1 23 23 23 ...
## $ WtavM : num 23.4 23.4 23.3 23.3 23.3 ...
## $ WtsdM : num 0.244 0.244 0.229 0.229 0.229 ...
## $ WtmaxR : num 24.1 24.1 24 24 24 ...
## $ WtminR : num 22.4 22.4 22.3 22.3 22.3 ...
## $ WtavR : num 23.1 23.1 23 23 23 ...
## $ WtsdR : num 0.424 0.424 0.417 0.417 0.417 ...
## $ WtrangeM : num 0.994 0.994 0.908 0.908 0.908 ...
## $ WtrangeR : num 1.76 1.76 1.75 1.75 1.75 ...
## $ StmaxM : num 32.3 32.3 32.3 32.3 32.3 ...
## $ StminM : num 31.7 31.7 31.8 31.8 31.8 ...
## $ StavM : num 32.2 32.2 32.1 32.1 32.1 ...
## $ StsdM : num 0.159 0.159 0.137 0.137 0.137 ...
## $ StmaxR : num 33.4 33.4 33.3 33.3 33.3 ...
## $ StminR : num 31 31 31.6 31.6 31.6 ...
## $ StavR : num 32.3 32.3 32.3 32.3 32.3 ...
## $ StsdR : num 0.545 0.545 0.432 0.432 0.432 ...
## $ StrangeM : num 0.634 0.634 0.534 0.534 0.534 ...
## $ StrangeR : num 2.39 2.39 1.71 1.71 1.71 ...
## $ winRSquared : num 0.874 0.874 0.865 0.865 0.865 ...
## $ winIntercept: num 21 21 20.9 20.9 20.9 ...
## $ winSlope : num 0.0959 0.0959 0.0938 0.0938 0.0938 ...
## $ ccoldday : int 42 42 42 42 42 42 42 42 42 42 ...
## $ sumRSquared : num 0.717 0.717 0.759 0.759 0.759 ...
## $ sumIntercept: num 32.5 32.5 32.6 32.6 32.6 ...
## $ sumSlope : num -0.0338 -0.0338 -0.0371 -0.0371 -0.0371 ...
## $ warmday : int 301 301 301 301 301 301 290 290 290 290 ...
Degree heating weeks is a measure of the cumulative effect of temperature (originally heat stress, but optionally cold stress) developed by NOAAs coral reef watch, with details here.
To calculate these metrics we will use the package dbcaDHW. DHW is calculated by comparing temperatures to a climatological threshold - The Maximum Monthly Mean (MMM) for each Site, which is calculated by averaging the annual max monthly means over a baseline climatology period (1985–1995).
# Read and prep data
temp <- read.csv("rawdata/ctemp_sst_mod.csv")
str(temp)
## 'data.frame': 12370 obs. of 11 variables:
## $ date : chr "1985JAN01" "1985JAN02" "1985JAN03" "1985JAN04" ...
## $ date1 : chr "1/1/85" "2/1/85" "3/1/85" "4/1/85" ...
## $ yearday : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Year : int 1985 1985 1985 1985 1985 1985 1985 1985 1985 1985 ...
## $ Delma : num 22.7 22.6 22.5 22.2 22.1 ...
## $ RasGhanada: num 23.3 23.1 23 22.8 22.8 ...
## $ Saadiyat : num 23.3 23.1 23 22.7 22.7 ...
## $ AlAqah : num 24 23.9 23.9 23.9 23.9 ...
## $ Dibba : num 24 23.9 23.9 23.8 23.8 ...
## $ month : chr "JAN" "JAN" "JAN" "JAN" ...
## $ season : chr "winter" "winter" "winter" "winter" ...
sites <- names(temp)[5:9]
g_num <- length(sites) + 1
temp1 <- temp[,c(2,5:9)]
temp1$date1 <- as.Date(temp1$date1, '%d/%m/%y')
## Note this is MODIFIED FROM ORIGINAL DHW FORMULA
## create a mean of the maximum monthly mean each year at each site from "past period"
mmtdata <- temp1 %>%
tidyr::gather("Site", "sst", 2:all_of(g_num)) %>%
#dplyr::filter(date1 < "2013-12-31") %>% #define past
dplyr::filter(date1 >= "1985-01-01" & date1 <= "1995-12-31") %>% #define past
dplyr::mutate(month = month(date1),
year = year(date1)) %>%
dplyr::group_by(Site, year, month) %>%
dplyr::summarise(avg = mean(sst)) %>% #obtain mthly avgs
dplyr::group_by(Site, year) %>%
dplyr::summarise(mmmt = max(avg)) %>% #obtain max mthly avg each year over this period
dplyr::group_by(Site) %>%
dplyr::summarise(mmmt = mean(mmmt)) #average this value as threshold for DHW
# The DHW calculation is assisted by creating a function to do the work of a rolling 12 week “sum”.
# A present period also needs to be defined.
## function for rolling cusum of 12 weeks - 84 days)
dhw_calc_12 <- tibbletime::rollify(sum, window = 84)
## observation period, hotspot & dhw calc
dhw <- temp1 %>%
tidyr::gather("Site", "sst", 2:all_of(g_num)) %>%
dplyr::filter(date1 >= "2000-01-01" & date1 <= "2013-12-31") %>% #define present
dplyr::mutate(month = month(date1)) %>%
dplyr::full_join(mmtdata, by = "Site") %>%
dplyr::mutate(hspt = sst - mmmt,
hsptm = ifelse(hspt >= 1, hspt, 0),
dhw = dhw_calc_12(hsptm)*(1/7)) %>% #hotspot based on greater than 1 degree as per coral watch
dplyr::select(-month, -hsptm)
dhwm <- dhw %>%
dplyr::mutate(year = year(date1)) %>%
dplyr::group_by(Site, year) %>%
dplyr::summarise(cummu_dhwm = sum(dhw, na.rm = TRUE))
The dhw object contains daily DWH values for each site. Because we need a summary of dhw by year, we can simply sum all DHW within each year as a yearly metric of total heat stress that year.
dhwm <- dhw %>%
dplyr::mutate(year = year(date1)) %>%
dplyr::group_by(Site, year) %>%
dplyr::summarise(cummu_dhwm = sum(dhw, na.rm = TRUE))
Here is what the DHW data look like for each site over the coral growth period at each site with NOAA bleaching alert levels indicated:
# Converting Site to factor for ordering
dhw <- dhw %>%
mutate(Site = as.factor(Site))
# Plot: DHW time series with alert thresholds
ggplot(dhw, aes(x = date1, y = dhw)) +
geom_line(color = "darkorange") +
geom_hline(yintercept = 4, linetype = "dashed", color = "blue") +
geom_hline(yintercept = 8, linetype = "dashed", color = "red") +
facet_wrap(~ Site, scales = "free_y", ncol = 1) +
labs(
title = "Degree Heating Weeks (DHW) 2000–2013",
subtitle = "Dashed lines: Bleaching Alert Levels (4°C-weeks = Blue, 8°C-weeks = Red)",
x = "Date",
y = "DHW (°C-weeks)"
) +
theme_minimal()
There is not a lot of variation in dhw, with most daily values throughout most years being zero. This is because for DWH at least, the reference period defining the threshold had generally greater temperatures compared to the period between 2000 and 2013.
The plot below shows the mean monthly temperature (MMT) within each year at each site. The red dashed line indicates the the maximum MMT and the blue dashed line the minimum MMT across all years at that site. The filled lines are the lowess fits illustrating the temporal trend in the data for the maximum MMT (red) and minimum MMT (blue). There are few dhw values > 0 in our dataset because there is a general trend that years (2000-2013) experienced lower maxmimum MMT compared to the reference period (1985-1995).
This plot also illustrates that in winter there is not as clear a temporal signal in minimum MMT, but there is arguably a small downtrend at a couple of sites.
Lets also calculate degree cooling weeks (DCW), where like DHW, the accumulated cold stress over a rolling 12-week (84-day) window is calculated, where daily sea surface temperature (SST) is at least 1°C below a site-specific cold threshold, and the excess cold is summed and expressed in °C-weeks
mmtdata_cool <- temp1 %>%
tidyr::gather("Site", "sst", 2:all_of(g_num)) %>%
#dplyr::filter(date1 < "2013-12-31") %>% #define past
dplyr::filter(date1 >= "1985-01-01" & date1 <= "1995-12-31") %>% #define past
dplyr::mutate(month = month(date1),
year = year(date1)) %>%
dplyr::group_by(Site, year, month) %>%
dplyr::summarise(avg = mean(sst)) %>% #obtain mthly avgs
dplyr::group_by(Site, year) %>%
dplyr::summarise(mmmt = min(avg)) %>% #obtain max mthly avg each year over this period
dplyr::group_by(Site) %>%
dplyr::summarise(mmmt = mean(mmmt)) #average this value as threshold for DHW
dhw_calc_12 <- tibbletime::rollify(sum, window = 84)
## observation period, hotspot & dhw calc
dhw_cool <- temp1 %>%
tidyr::gather("Site", "sst", 2:all_of(g_num)) %>%
dplyr::filter(date1 >= "2000-01-01" & date1 <= "2013-12-31") %>% #define present
dplyr::mutate(month = month(date1)) %>%
dplyr::full_join(mmtdata_cool, by = "Site") %>%
dplyr::mutate(hspt = sst - mmmt,
hsptm = ifelse(hspt <= -1, hspt, 0),
dcw = dhw_calc_12(hsptm)*(1/7)) %>% #hotspot based on less than 1 degree
dplyr::select(-month, -hsptm)
## again, because we need a summary of dhw by year, we can simply sum all DHW within each year as a yearly metric of total heat stress that year.
dhwm_cool <- dhw_cool %>%
dplyr::mutate(year = year(date1)) %>%
dplyr::group_by(Site, year) %>%
dplyr::summarise(cummu_dcwm = sum(dcw, na.rm = TRUE))
And a plot of DCW values over time: